Momentum distribution for a one-dimensional trapped gas of hard-core bosons 
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I. INTRODUCTION 

Recent advances in atom waveguide technology [[T UTUll and 
the realization of Bose-Einstein condensates in optical and 
magnetic traps of variable aspect ratio fTTl fl2l have spurred 
interest in the properties of degenerate quantum gases in lower 
dimensions. In particular, the Tonks gas, in which strong 
transverse confinement and low temperature and density al- 
low the gas to be modeled as a one-dimensional (ID) system 
of point particles with hard-core interactions lfT3l FBI , is of 
considerable theoretical interest due to the fact that it defies a 
mean-field description, but is on the other hand exactly solu- 
ble via the Fermi -Bose mapping |[T5l[T6ll . Although the exact 
many-body wave function can be written in a compact form 
using the mapping theorem, obtaining information about im- 
portant observables has proven to be a difficult task. One such 
quantity that bears the signature of the Tonks gas is the mo- 
mentum distribution, which has a sharp peak at zero momen- 
tum IT3l . in contrast to the Fermi sea for the corresponding 
ID gas of fermions. In a mathematical tour de force, Lenard 
ifTTl obtained upper bounds on the momentum distribution for 
a homogeneous Tonks gas, and some elaborations of that work 
followed ifTHIl . In a previous paper [19| we obtained numer- 
ical results for the momentum distribution of a harmonically 
trapped Tonks gas for up to A" = 10 particles, and Minguzzi 
et al. [20| developed an analytic approximation for the high- 
momentum tail of the momentum distribution of a trapped 
gas. Cazalilla ETI obtained an analytic approximation for the 
momentum distribution of Tonks gas confined in a box using 
a description of the system as a Luttinger liquid. 

Our previous calculations of a trapped Tonks gas were per- 
formed using a Monte Carlo (MC) integration of the many- 
body wave function to obtain the single-particle reduced den- 
sity matrix from which the momentum distribution was ob- 
tained via Fourier transformation. Although the data thus gen- 
erated were useful, an improved method is desirable because 
of the limited accuracy of MC integration. Even these MC 
data were limited to A^ = 10 particles, which required weeks 
of computer time. It seems clear that the way forward is to 
have high-precision numerical data available for testing the 
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validity of approximations. In this paper we start from the 
A^-particle ground state wave function for a one-dimensional 
condensate of hard-core bosons in a harmonic trap and de- 
velop an algorithm to compute the reduced single-particle 
density matrix and momentum distribution. The key advan- 
tage of this approach is that, although we are limited to only 
N = 8 particles at present by computer resources, these data 
are of high precision, thus providing a testing ground for ana- 
lytic approximations. 

In Sec. [Il] we give a precise definition of the system and 
find the ground state wave function using the Fermi-Bose 
mapping theorem. In Sec. Ill we write the single-particle 



reduced density matrix p\ and develop a method for its nu- 
merical solution. In Sec. [IV] we use these results for p\ to 
evaluate the momentum distribution and compare the results 
to a recent approximation for the high-momentum tail. 



H. GROUND STATE WAVE FUNCTION 
The Hamiltonian of A" bosons in a ID harmonic trap is 

h 2 d 2 l 



N 



2m dx 2 
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We assume that the two-body interaction potential consists 
only of a hard-core of ID diameter a. This is conve- 
niently treated as a constraint on allowed wave functions 
ip(xi, . . . , xn) such that 



■0 = if \xj - Xk\ < a, l<j<k<N, 



(2) 



rather than as an infinite interaction potential. It follows from 
the Fermi-Bose mapping theorem [ il31 [161 l22l that the exact 
A^-boson ground state ipso of the Hamiltonian (1) with the 
constraint (2) is 



^Bo(Xi,...,X N ) = \'ip F o(Xl,...,X N )\, 



(3) 



where 4>fo is the ground state of a fictitious system of A" spin- 
less fermions with the same Hamiltonian ([1} and constraint. 
At low densities it is sufficient lfl3l [Pfl to consider the case 
of impenetrable point particles, the zero-range limit a — > of 
Eq. (2). Since wave functions of "spinless fermions" are anti- 
symmetric under coordinate exchanges, their wave functions 
vanish automatically whenever any ij — Xk, the constraint 



2 



has no effect, and the corresponding fermionic ground state is 
the ground state of the ideal gas of fermions, a Slater determi- 
nant of the lowest N single-particle eigenfunctions </>„ of the 
harmonic oscillator (HO) 



I (N-l,N) 

^fo(xi, ■ ■ ■ ,X N ) = — r == det i n {Xj). 

VNl («j)=(o,i) 



(4) 



The HO orbitals are 

4>n{x) 



e- x2 ' 2 H n {x/ Xosc ), (5) 



with H n {x) the Hermite polynomials and x osc — ^Jhjmuj 
the ground state width of the harmonic trap for a single atom. 
For convenience, we introduce the dimensionless coordinates 
%i = £i/x osc , and the corresponding ground state wave func- 
tion ipso- As we have shown in previous work |fl9l , substi- 
tution of Eq. <|3j into Eq. (|4]) and some matrix manipulations 
ll23l lead to a simple but exact expression of the Bijl-Jastrow 
pair product form for the TV-boson ground state: 



^bo{xi, ■ . . ,Xn) = C 



N 



' N 

.1=1 



n \ x k-xj\, 

l<j<k<N 

(6) 



with normalization constant 



C N = 2 N ^ N ^/ A 



N-l 



n=0 



-1/2 



(7) 



III. SINGLE-PARTICLE DENSITY MATRIX 

A. Analytic formula 

The reduced single-particle density matrix with normaliza- 
tion J pi(x, x)dx — N for the ground state given by Eq. (|6jl 



pi(x,x f ) = N J ip B0 (x,x 2 ,...,x N ) 

xip B0 (x', x 2 , ■ ■ ■ ,x N )dx 2 ■ ■ -dx N 
= M N e- x2 l 2 e- x ' 2 ' 2 I{x,x'), (8) 

where the integration is from — oo to oo for each coordinate 
unless otherwise stated, and 



A/V = iV2 7V ^- 1 )/ 2 7r- JV / 2 



' N 

n=0 



(9) 



and we have defined 



r N 

I(x,x) = I TT e~ Xi \xi — x\\xi — x\ 



i=2 

x Y\ (xk - Xj) 2 dx 2 ■ ■ ■ dxN- (10) 

2 <j<k<N 



In the following subsection, we develop a method for analyz- 
ing I(x, x'). We will see that, when N is small enough (say 
2 or 3), the exact expression is manageable, but that we must 
turn to numerical methods for larger N. 



B. Numerical approach 



The multidimensional integral (lOi can be expressed in 



terms of polynomials, Gaussians, and error functions. But, 
even for relatively small N, the number of terms in such an 
expression is too large to be useful when written. We pre- 
viously evaluated the reduced single-particle density matrix 
using MC methods lfl9l . but with limited numerical accuracy. 
Here we present a seminumerical approach in which we repre- 
sent the integral in terms of incomplete gamma functions and 
evaluate the result to machine precision for particular values 
of x and x' . 

We reduce the integral I(x, x') to incomplete gamma func- 
tions in the following way. Consider the case x < x' . We first 
integrate over x 2 , writing 



I(x,x') 



f N 

/ dx$ ' ■ ■ dx j\[ I 2 J^J g x * \xi x'l 
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x n (Zfc-Zj) 2 ; 
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where 



h = 



I 2 (x,x',x 3 , . . .,X N ) 

I - + \P 2 dx 2 , x' <x (12) 

I J — oo J x J x' 



and 



P 2 = P 2 (x,x',x 2 ,...,x N ) 

= e - x >{x 2 -x){x 2 -x') H (x 2 -x k f. (13) 

2<fc<7V 

The integrand P 2 is an analytic function of x, x', and x 2 ( in 
the sense that derivatives of all orders in these variables exist), 
so that the integral over each of the three intervals is analytic 
in these variables. Furthermore, we can evaluate the integral 
over x 2 in Eq. ( |T2] > easily because the integrand is a Gaussian 
multiplied by a polynomial. We next integrate over 2:3, getting 

N 



I{x, x') = / dx4 ■ ■ ■ dx^h Y\ 



G 1 I X 1 ^11 X % X 

x n ( x k~xj) 2 , 

4<j<k<N 



(14) 



where 



h = h(x,x',X4,...,x N ) 

( px px 1 pOO j 

= I - + \hP3dx 3 , x 1 <x (15) 

I J —00 J x J x' 
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FIG. 1: Gray scale plots of the dimensionless reduced density matrix x osc pi(x, x') as a function of the dimensionless coordinates x and x' , 
for (a) N = 2, (b) N = 6, and (c) N = 8. 



and 



Pz = P3(X,X',X3,...,X N ) 

= e~ x 3{x 3 - x)(x 3 - x') Yl (x 3 -x k ) 2 . (16) 

3<k<N 

It is important to note that, because I2 is a polynomial in 
x 3 , . . . ,Xn, the integral ^3 in Eq. ( 15 1 can be solved by the 
same technique used to solve I2 in Eq. (12 1. We continue 



this procedure, defining I4, etc., until all of the integrals are 
finished. At each stage, one has a more complicated polyno- 
mial in the remaining independent variables in the integrand. 
Before continuing, we note that there is, of course, nothing 
essentially new when we choose x' < x. For instance, for I2, 
we have 




P2 dx2, x' < x. 



(17) 



If we examine each of the integrals I2 , 13, . . ■ above, we see 
that all of the integrals to be computed can be reduced to terms 
proportional to integrals of the form 



to compute ^3, and so on, until all powers of x 7 ^ have been 
replaced by a n (x,x'). This procedure can be simplified by 
performing all of these substitutions at once. To this end, con- 
sider 



JV 



P = Y\_(xi ~ x)(xi 



i=2 



x') n (x 3 

2<j<k<N 



Xk) 2 



(20) 



To compute the integral I(x,x') given in Eq. (10 1, we ex- 
pand Eq. (20 1 and substitute a n (x, x') for each occurrence of 
x^, for any m. In addition, for each value of m for which a 
term is independent of x m , a factor of ao(x,x') must be in- 
cluded. The result is I(x, x') expressed as a polynomial in the 
a n (x,x') for approximately 2N values of n. Rather than print 
the results, we store a table of the coefficients of the powers 
of a n (x,x') on a computer. Then, a table of the a n (x,x') 
is computed for a particular pair x, x', and I(x, x') is com- 
puted using this table together with the table of coefficients. 
Evaluating the a n (x,x') using numerical integration is rela- 
tively inefficient. Instead, we evaluate them numerically using 
well-known, efficient routines to compute incomplete gamma 
functions. We define an indefinite integral 



x px />oo 



&n (Xj X ) 



— 00 J x J x 



and 



i(x,x') = 




+ / J> e~ u u n du, x < x' 

(18) 



du x > x 1 



(19) 

We now present an algorithm for expressing the integrals 
I2, ■ ■ ■ in terms of the a n (x,x') defined in Eqs. (I81 and 
( |19) , To achieve this for I 2 , we take P 2 , drop the factor of 
exp(— X2), expand the remaining factors as a polynomial in 
X2, and replace each occurrence of x r 2 L with a n (x,x'). The 
result is I2 expressed as a polynomial in x, x',x 3) . . . with co- 
efficients involving the a n (x, x'). This expression for I 2 is 



substituted into Eq. ( 15 1 and the same procedure is then used 



F n (x) = I e~ x 'x n 



(21) 



Then using definitions ( fl~8] > and ( [T9| we have 



a n (x,x') = F„(oo) - F„(-oo) 

+ sga(x' -x)2[F n (x)-F n (x')}. (22) 

a n (x, x') is continuous, but has a cusp at x = x'. For numer- 
ical computation, we use 




if n is odd , 

, . . . _(n-l)!! (23) 
r I J = 2 „ /2 otherwise, 
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and 



2[F n (x) - F n (x')] = [sgn(x)] n r ( — ,x 2 



sgn.:,-';."l ( ' .;".,•'-•) : • 



C. Examples 

In this section we carry out in detail the algorithm given 
above for two particles and give the result for three particles. 
We first choose N = 2 and tabulate the required values of F n 
and a n (x, x'). Integrating Eq. (21 1 by parts, we write 



1 1 _ 2 

Fo(x) — -\Zirerfx xe x 

ZK ; 4 2 



Using Eq. d22b we then have 



a (x,x') 
ai(x, x 1 ) 

a 2 {x, x') 



\fn + sgn(a/ — x)\/tt (erfa — erfa') , 
-sgn(x' - x) (e~ x2 - e~ x 



+ sgn(x' — x) 



(erfa — erfa') 



+x e — xe 



Fo(x) 
F 1 (x) 



-erfa, 



Applying the algorithm outlined in the previous section we 
find for two particles 



3-2 X I2 

pi{x,x') = M 2 e 2 2 



x 2 x' 2 

e~ X2 \(x2 — x)(x 2 — x')\ dx2 — N 2 e^^ »" [0:2(2;, x 1 ) — (x + x')ai(x, x') + xx'cto{x 7 x')] 



Af 2 e 2 2 < \pK 



2 +xx 



sgn(a;' — x) 



- + xx' I \Pk ( erfa — erfa' 



For N = 3, one can easily compute the expression for 
I(x, x') in terms of a n (x, x') by hand, with the result 

x 2 x' 2 

pi{x,x) = 2A/" 3 e _ T 2- { _ a \ + a 2 a 4 

+ p 2 {—a\ + a\az) + t(a 2 — 2a\a 3 + aoa^) 
+ t 2 {-a\ + a a 2 ) 

+ p[a 2 a^ — aiot4 + t{a\a 2 — aoo^)] }, 

where p = x + x', t = xx', and the explicit dependence of 
a n on x and x' is omitted for clarity. For larger values of N, 
I(x,x') rapidly becomes more difficult to compute by hand. 
In the next section we present the results of carrying out the 
algorithm detailed above on a computer. 



D. Numerical results 

We have evaluated the above integrals numerically for N — 
2-8. 

Figure [T] shows a gray scale plot of the dimensionless re- 
duced single-particle density matrix x osc pi(x, x') versus the 
normalized coordinates x and x' for (a) N = 2, (b) N = 6, 
and (c) N — 8. We verified that along the diagonal pi (x,x) = 
p(x) reproduced the single-particle density l24l . 



IV. MOMENTUM DISTRIBUTION 

In terms of the boson annihilation and creation operators 
in position representation (quantized Bose field operators) the 
one-particle reduced density matrix is 

pi(x,x') = (*B0$V)tol*B0>. (25) 

The momentum distribution function n(k), normalized to 
JZo n ( k ) dk = N ' is "( fc ) = (^Bo\aHk)a{k)\^ BO ) where 
a(k) is the annihilation operator for a boson with momentum 
hk. Then 

/OO /"OO 
dx / dx' Pl (x,x')e~ lk{x - x,) . (26) 
-co J — OO 

The spectral representation of the density matrix then leads to 
n(k) = J^j \k i 3 (k) | 2 where the [ij are Fourier transforms of 
the natural orbitals: p,j(k) = (2ir)~ 1 ^ 2 (j) n (x)e~ lkx dx. 
Figure|2]shows the numerically calculated dimensionless mo- 
mentum distribution k osc n(K) versus normalized momentum 
k = fc/fc osc , with fc osc = 2ir/x osc , for (a) N = 2, (b) N = 6, 
and (c) N = 8. We typically evaluated p 1 (x,x') to machine 
precision for smaller values of k, with precision decreasing to 
a part in 10 -6 for the largest values of k. The key features are 
that the momentum distribution maintains the peaked struc- 
ture reminiscent of the spatially uniform case lfl3l l25l for the 
ID HO, and that the peak becomes sharper with increasing 
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FIG. 2: Dimensionless momentum distribution fc osc n(re) versus nor- 
malized momentum k — k/k osc for N = 2, N — 6, and N — 8. 
Note the peaks becoming sharper with increasing N. 
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atom number N. This is to be expected since as the number 
of atoms increases the many-body repulsion causes the sys- 
tem to become more spatially uniform within the trap interior. 
Minguzzi et al. [20 | determined that the momentum distribu- 
tion for finite N decays according to 



lim k 4 n(k) = An, 



(27) 



where depends only on the number of particles. In partic- 
ular, for N — 2, they found 



-4, 



(Smu)- 3/a . 



(28) 



K 



Figure [3] shows our numerical results for N = 2 approach- 
ing this asymptotic form. The dashed line in Fig. [3] shows 
kg SC A2/ 'k 4 versus k = k/k osc , and we see that this approxi- 
mation agrees with our numerical results for N = 2 for high 
momenta. Furthermore, inspection of the numerical results 
for other values of N shows a 1/fc 4 dependence in the high- 
momentum tail. 

V. SUMMARY AND CONCLUSIONS 

In summary, we have developed a method for obtain- 
ing high-accuracy results for the momentum distributions of 
trapped Tonks gases, and presented results for up to eight par- 
ticles. Our results agree reasonably with the high-momentum 
approximation n(p) oc l/p 4 obtained by Minguzzi et al. 
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FIG. 3: Dimensionless momentum distribution fc osc n(re) versus nor- 
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